# With Village Fixed Effects ---------------------------------------------------

# Remove observations that contain missing values for electricity supply dimensions
pool_elec <- pool %>% dplyr::select(elec_respvars, village, hh, age, hhsize, bplaay,
                                    bplaay, wave, gender, religion, education, lnexp,
                                    weights) %>% 
        na.omit()

# Duration Day (Total hours per day)
model.dd.vfe.duration_day.bplaay = fixest::feols(duration_day ~ bplaay*wave + age + gender +
                                                 religion + education + hhsize +
                                                 lnexp  | village, 
                                               data = pool_elec, weights = pool_elec$weights)

# Duration Night (Hours in the evening)
model.dd.vfe.duration_night.bplaay = fixest::feols(duration_night ~ bplaay*wave + age + gender +
                                                   religion + education + hhsize +
                                                   lnexp  | village, 
                                                 data = pool_elec, weights = pool_elec$weights)

# Reliability (Outage days per month)
model.dd.vfe.reliability.bplaay = fixest::feols(reliability ~ bplaay*wave + age + gender +
                                                religion + education + hhsize +
                                                lnexp  | village, 
                                              data = pool_elec, weights = pool_elec$weights)

# Voltage Fluctuation (Voltage fluctuation days per month)
model.dd.vfe.quality_fluc.bplaay = fixest::feols(quality_fluc ~ bplaay*wave + age + gender +
                                                 religion + education + hhsize +
                                                 lnexp  | village, 
                                               data = pool_elec, weights = pool_elec$weights)

# Low Voltage (Low Voltage days per month)
model.dd.vfe.quality_low.bplaay = fixest::feols(quality_low ~ bplaay*wave + age + gender +
                                                religion + education + hhsize +
                                                lnexp  | village, 
                                              data = pool_elec, weights = pool_elec$weights)

# With Household Fixed Effects -------------------------------------------------

# Duration Day (Total hours per day)
model.dd.hfe.duration_day.bplaay = fixest::feols(duration_day ~ bplaay*wave + hhsize +
                                                 lnexp  | hh, 
                                               data = pool_elec, weights = pool_elec$weights)

# Duration Night (Hours in the evening)
model.dd.hfe.duration_night.bplaay = fixest::feols(duration_night ~ bplaay*wave + hhsize +
                                                   lnexp  | hh, 
                                                 data = pool_elec, weights = pool_elec$weights)

# Reliability (Outage days per month)
model.dd.hfe.reliability.bplaay = fixest::feols(reliability ~ bplaay*wave + hhsize +
                                                lnexp  | hh, 
                                              data = pool_elec, weights = pool_elec$weights)

# Voltage Fluctuation (Voltage fluctuation days per month)
model.dd.hfe.quality_fluc.bplaay = fixest::feols(quality_fluc ~ bplaay*wave + hhsize +
                                                 lnexp  | hh, 
                                               data = pool_elec, weights = pool_elec$weights)

# Low Voltage (Low Voltage days per month)
model.dd.hfe.quality_low.bplaay = fixest::feols(quality_low ~ bplaay*wave + hhsize +
                                                lnexp  | hh, 
                                              data = pool_elec, weights = pool_elec$weights)

## Model and coefficient lists for village and household fixed effects combined ##
models.vfe.dd_elec_bplaay <- list(model.dd.vfe.duration_day.bplaay,model.dd.vfe.duration_night.bplaay,
                                model.dd.vfe.reliability.bplaay, model.dd.vfe.quality_fluc.bplaay,
                                model.dd.vfe.quality_low.bplaay)

models.hfe.dd_elec_bplaay <- list(model.dd.hfe.duration_day.bplaay, 
                                model.dd.hfe.duration_night.bplaay, model.dd.hfe.reliability.bplaay, 
                                model.dd.hfe.quality_fluc.bplaay, model.dd.hfe.quality_low.bplaay)

# Generate latex tables
esttex(models.vfe.dd_elec_bplaay, 
       se = "cluster", 
       digits = 3,
       dict = c(bplaay = "BPL/AAY", wave = "2018", duration_day = "Daily Supply Hours",
                duration_night = "Evening Supply Hours", reliability = "Monthly Outage Days",
                quality_fluc = "Fluctuating Voltage Days", quality_low = "Low Voltage Days",
                village = "Village", hh = "Household"),
       keep = c(":", "2018"),
       fitstat = c("r2"),
       fixef_sizes = T,
       float = F,
       replace = T,
       file = c("./Manuscript/Tables/dd_elec_vfe_bplaay.tex"))

esttex(models.hfe.dd_elec_bplaay, 
       se = "cluster", 
       digits = 3,
       dict = c(bplaay = "BPL/AAY", wave = "2018", duration_day = "Daily Supply Hours",
                duration_night = "Evening Supply Hours", reliability = "Monthly Outage Days",
                quality_fluc = "Fluctuating Voltage Days", quality_low = "Low Voltage Days",
                village = "Village", hh = "Household"),
       keep = c(":", "2018"),
       fitstat = c("r2"),
       fixef_sizes = T,
       float = F,
       replace = T,
       file = c("./Manuscript/Tables/dd_elec_hfe_bplaay.tex"))
